function U = LTVD(Uin, Ori, lambda,method)
% LTVD - TVD格式的迭代步
%   INPUTS:
%       Uin: 前一步的值。1*n的cell数组
%       Ori: Orientation表示这个tvd的方向x,或者y,用于确定边界条件。
%           1：x向，且右端出流，台阶上面。
%           2：x向，且右端在台阶下。
%           3：y向。
%       lambda: lambda = tau/h;
%       method:方法选项
%           1：ULT1(default)
%           2：ULT1C
%    CALL:
%    LTVD(Uin,boundary,lambda,method)
    n= length(Uin);m=4; 
    U_hr=zeros(m,n+3);                  % U_hr(j+2)=U_{j+1/2}
    U_go=zeros(m,n+4);                  % hr is short for half and right
    for i = 1:n
        U_go(:,i+2)=Uin{i};             % U_go(j) = U_{j-2}
    end
    if Ori == 1 || Ori ==2              % 入流
                                        %U_go(:,1)=[1.4;3;0;6.55];%
        U_go(:,1)=Variables2U(1.4,3,0,1);
        U_go(:,2)=U_go(:,1);
    else                                % y向反射，v取反向。
        U_go(:,2)=U_go(:,3);
        U_go(3,2)=-U_go(3,2);           % v=-v
        U_go(:,1)=U_go(:,4);
        U_go(3,1)=-U_go(3,1);
    end

    if Ori == 1                         % 出流边界。偏导为0
        U_go(:,n+3)= U_go(:,n+2);
        U_go(:,n+4)= U_go(:,n+2);
    else      
        U_go(:,n+3)= U_go(:,n+2);
        U_go(:,n+4)= U_go(:,n+1);
        if Ori == 2                     % x向反射边界
            U_go(2,n+3)= -U_go(2,n+2);  % u = -u
            U_go(2,n+4)= -U_go(2,n+1);  % u = -u
        else                            % y向反射边界
            U_go(3,n+3)= -U_go(3,n+2);  % v = -v
            U_go(3,n+4)= -U_go(3,n+1);  % v = -v
        end
    end
    myplot(U_go,'U ghost');
    
    U_hr = 0.5*(U_go(:,1:(n+3))+U_go(:,2:(n+4)));
    myplot(U_hr,'U j+1/2');
    

    R = cell(1,n+3);                    % R(j+2)=R_{j+1/2}
    L = cell(1,n+3);                    % L(j+2)=L_{j+1/2}
    nu = zeros(m,n+3);                  % \nu
    for i = 1:(n+3)
        v = U_hr(:,i);
        %[a,L{i},R{i}]=Eigen(v,Ori);
        %nu(:,i)=lambda*a;
        if Ori == 1||Ori == 2           % x向
            [a,R{i}]=Rx(v);             % R_{j+1/2}=Rx(U_{j+1/2})
            nu(:,i)=lambda*a;
            L{i}=Lx(v);                 % L_{j+1/2}=Lx(U_{j+1/2})
        else                            % y向
            [a,R{i}]=Ry(v);             % R_{j+1/2}=Ry(U_{j+1/2})
            nu(:,i)=lambda*a;           % 
            L{i}=Ly(v);                 % L_{j+1/2}=Ly(U_{j+1/2})
        end
    end
    alpha = zeros(m,n+3);               % alpha(j+2)=alpha_{j+1/2}
    for i = 1:(n+3)
        v = U_go(:,i+1)-U_go(:,i);
        alpha(:,i)=L{i}*v;              % alpha_{j+1/2}=L_{j+1/2}*(U_{j+1}-U_{j})
    end
    myplot(alpha,'\alpha')
    
    if method == 2
        theta = zeros(m,n+2);               % theta(j)=theta_{j-1},
        for i=1:(n+2)                       % 
            for k = 1:m
                if alpha(k,i+1)~=0 ||alpha(k,i)~=0
                    theta(k,i)=abs(alpha(k,i+1)-alpha(k,i));
                    theta(k,i)=theta(k,i)/(abs(alpha(k,i+1))+abs(alpha(k,i))); 
                else
                    theta(k,i)=0;
                end            
                % theta_{j-1}=abs(alpha_{j-1/2}-alpha_{j-3/2})/(|alpha_{j-1/2}|+|alpha_{j-3/2}|)
            end
        end
        htheta = zeros(m,n+1);              % htheta(j+1)=\hat{theta}(j+1/2)
        for i=1:(n+1)
            for k = 1:m
                htheta(k,i)=max(theta(k,i),theta(k,i+1));
            end
        end
    end
    
    tg = zeros(m,n+3);                  % tg(j+2)=\tilde{g}_{j+1/2}
    for i =1:(n+3)
        for k = 1:m
            tg(k,i)=( Q( nu(k,i) ) - nu(k,i)^2 )*alpha(k,i)/2;
        end
    end
    myplot(tg,'\tilde{g}')
    g = zeros(m,n+2);                   % g(j+1)=g_{j}
    for i=1:(n+2)
        for k=1:m
            s = sign(tg(k,i+1));
            g(k,i)=s*max(0,min( abs(tg(k,i+1)), s*tg(k,i) ));
        end
    end
    myplot(g,'g')
    gamma = zeros(m,n+1);               % gamma(j+1)=gamma_{j+1/2}
    for i = 1:(n+1)
        for k = 1:m
            if alpha(k,i+1) ==0
                gamma(k,i)=0;
            else            
                gamma(k,i)=(g(k,i+1)-g(k,i))/alpha(k,i+1);
            end
        end
    end
    myplot(gamma,'\gamma')
    beta = zeros(m,n+1);                % beta(j+1)=beta_{j+1/2}
    for i = 1:(n+1)
        for k = 1:m
            if method == 2
                s = (1+ 2*htheta(k,i) );
            else
                s=1;
            end
            beta(k,i)=Q( nu(k,i+1) + s*gamma(k,i) )*alpha(k,i+1)-s*(g(k,i)+g(k,i+1));
        end
    end
    myplot(beta,'\beta')
    f = zeros(m,n+2);
    if Ori == 1||Ori == 2               % x向
        for i=1:(n+2)
            f(:,i) =F(U_go(:,i+1));
        end
    else                                % y向
        for i=1:(n+2)
            f(:,i) =G(U_go(:,i+1));
        end
    end
    
    tf = zeros(m,n+1);
    for i = 1:(n+1)
        tf(:,i)=(f(:,i)+f(:,i+1));
        temp = zeros(m,1);
        for k = 1:m
            temp = temp+beta(k,i)*R{i+1}(:,k);
        end
        tf(:,i)=tf(:,i)-temp/lambda;
        tf(:,i)=tf(:,i)/2;
    end
    myplot(tf,'\tilde(f)');
    U = cell(1,n);
    for i = 1:n
        U{i}=Uin{i}-lambda*( tf(:,i+1) - tf(:,i) );
    end
end
function myplot(U,t)
% MYPLOT - debug
%   
%    hold on; %set(gca,'DataAspectRatio',[1 1 1]);set(gca,'Xtick',0:0.2:3);
             %    contourf(x,y,rho);
    return
    U
    plot(U(1,:))
    title(t);
    hold on
    plot(U(2,:))
    plot(U(3,:))
    plot(U(4,:))
    hold off;
    pause
end
